home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / SeqPups / appsrc / autoseq.src / CTraceFile.C < prev    next >
Text File  |  1996-07-05  |  45KB  |  1,406 lines
  1. //    ============================================================================
  2. //    CTraceFile.C                                                    80 columns
  3. //    Reece Hart    (reece@ibc.wustl.edu)                                tab=4 spaces
  4. //    Washington University School of Medicine, St. Louis, Missouri
  5. //    This source is hereby released to the public domain.  Bug reports, code
  6. //    contributions, and suggestions are appreciated (to email address above).
  7. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  8. //    Please see CTraceFile.H for a description of the CTraceFile class.
  9. //
  10. //    NOTES
  11. //    *    I make frequent use of the following construct:
  12. //            return (error = <error_t>);
  13. //        This causes the error flag of this instance of CTraceFile to be set
  14. //        to the appropriate error, and then returns this value.
  15. //        although CPeakList contains a lot that's irrelevant to CTraceFile.
  16. //    *    WriteSCF should scale the traces for writing to [0,255] range without
  17. //        affecting the trace permanently.
  18. //    ========================================|===================================
  19. #include    "CTraceFile.H"
  20. #include    <iostream.h>
  21. #include    <fstream.h>
  22. #include    <string.h>
  23. #include    <stdio.h>
  24. #include    <ctype.h>
  25. #include    <stdlib.h>
  26. #include    <time.h>
  27. #include    "CPeakList.H"
  28. #include    "RInclude.H"
  29. #include    "RInlines.H"
  30. #include    "DNA.H"
  31.  
  32. vrsn CTraceFile::version = "94.05.03";
  33.  
  34. const bool DEBUG = !TRUE;
  35.  
  36. //    ============================================================================
  37. //    CTraceFile (constructor)
  38. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  39. //    This constructor takes up to two arguments which specify filename and format
  40. //    for initialization.  If neither is provided, the instance is initialized
  41. //    with empty values.  If a filename is provided, then an attempt is made to
  42. //    read the file in the format specified by fmt, or in a format inferred by
  43. //    Read if fmt is not specified.  ie. The following are valid:
  44. //        CTraceFile    myTrace();            // empty instance
  45. //        CTraceFile    myTrace("myFile");    // attempt read myFile, guessing format
  46. //        CTraceFile    myTrace("myFile",SCF);    // attempt to read myFile in SCF fmt
  47. //    ========================================|===================================
  48. CTraceFile::CTraceFile(const char* fn, format_t fmt):
  49.     filename(NULL),
  50.     error(noError),
  51.     nativeFormat(unknown),
  52.     whichStrand(top),
  53.     numPoints(0),
  54.     min(0),
  55.     max(0),
  56.     mean(0),
  57.     leftCutoff(0),
  58.     rightCutoff(0),
  59.     numBases(0),
  60.     sequence(NULL),
  61.     basePositions(NULL),
  62.     comments(NULL)
  63.     {
  64.     for (uint traceIndex=0; traceIndex<NUM_BASES; traceIndex++)
  65.         trace[traceIndex]=NULL;
  66.  
  67.     if (fn)
  68.         error = Read(fn,fmt);
  69.     }
  70.  
  71. //    ============================================================================
  72. //    ~CTraceFile (destructor)
  73. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  74. //    Call Deallocate to clean up our dynamically-allocated space.
  75. //    Note that the lifetime of peaks (a CSequence) is the same as the lifetime
  76. //    of CTraceFile because it is automatically allocated when an instance of
  77. //    CTraceFile is allocated.  Therefore, it will be disposed of automatically
  78. //    its destructor will free all space allocated for the peak records.
  79. //    ========================================|===================================
  80. CTraceFile::~CTraceFile(void)
  81.     {
  82.     Deallocate();
  83.     }
  84.  
  85. //    ============================================================================
  86. //    Allocate
  87. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  88. //    Allocate space for the traces, bases, base_positions, and comment.  This is
  89. //    all-or-none: on success, we've allocated all (and returned noError); on
  90. //    failure we've allocated nothing (and returned memError);
  91. //    ========================================|===================================
  92. CTraceFile::error_t
  93. CTraceFile::Allocate(size_t pt_sz, size_t base_sz, size_t comment_sz)
  94.     {
  95.     // Allocate space for trace[]
  96.     for (uint traceIndex=0; traceIndex<NUM_BASES; traceIndex++)
  97.         if ( NULL == (trace[traceIndex]=new CTrace<trace_t>(pt_sz)) )
  98.             // allocate failed
  99.             {
  100.             Deallocate();
  101.             return (error = memError);
  102.             }
  103.         
  104.     // Allocate space for sequence and base_positions
  105.     if (    !(sequence = new char[base_sz]) ||
  106.             !(basePositions = new tracepos[base_sz]) ||
  107.             !(comments = new char[comment_sz+1]) )
  108.         {
  109.         Deallocate();
  110.         return (error = memError);
  111.         }
  112.     
  113.     numPoints = pt_sz;
  114.     numBases = base_sz;
  115.     rightCutoff = pt_sz-1;
  116.  
  117.     return (error = noError);
  118.     }
  119.  
  120. //    ============================================================================
  121. //    Deallocate
  122. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  123. //    Make Mom proud and put away all of our toys... ensure that all of our space
  124. //    is neatly reclaimed.  This method is provided separately so that in the case
  125. //    of errors, we can deallocate stuff in one call.  Note that this doesn't
  126. //    deallocate peaks; the reasoning is twofold:
  127. //    1) if we're ready to calculate the peaks, then we've already got our data
  128. //        allocated in its entirety. Thus, there's no reason to call this except
  129. //        to cause problems.  You wouldn't want to do that.
  130. //    2) It's probably not kosher to deallocate peaks... it's been allocated
  131. //        automatically and doing so might wreak havoc.
  132. //    ========================================|===================================
  133. void
  134. CTraceFile::Deallocate(void)
  135.     {
  136.     delete filename;
  137.     for (uint traceIndex=0; traceIndex<NUM_BASES; traceIndex++)
  138.         delete trace[traceIndex];
  139.     delete sequence;
  140.     delete basePositions;
  141.     delete comments;
  142.     }
  143.  
  144. //    ============================================================================
  145. //    CalculateStats
  146. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  147. //    Calculates the stats for individual traces and then calculates the overall
  148. //    stats for the combination of the four.  Theoretically, sum may overflow,
  149. //    although it's unlikely... if we have that many values, something's amiss.
  150. //    ========================================|===================================
  151. void
  152. CTraceFile::CalculateStats(void)
  153.     {
  154.     double        sum = 0;
  155.     
  156.     for (uint traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  157.         {
  158.         trace[traceIndex]->CalculateStats();
  159.  
  160.         if ((traceIndex == A) ||                // this is first pass ||
  161.             (trace[traceIndex]->Max() > max))    //    Max() for this trace > max
  162.             max = trace[traceIndex]->Max();        // set the new max
  163.  
  164.         if ((traceIndex == A) ||                // ditto for min values
  165.             (trace[traceIndex]->Min() < min))
  166.             min = trace[traceIndex]->Min();
  167.  
  168.         sum += trace[traceIndex]->Mean();        // cumulative value for means
  169.         }
  170.  
  171.     mean = (trace_t)sum/NUM_BASES;            // mean over 4 traces
  172.     }
  173.  
  174. //    ============================================================================
  175. //    Transform
  176. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  177. //    Applies a 4x4 transformation/orthogonalization matrix to the 4 traces,
  178. //    producing a new set of traces which replaces the existing set.
  179. //    Transformation matrics are expected to be 4x4 matrix of the form:
  180. //        M =            [,A]    [,C]    [,G]    [,T]
  181. //            [A,]    mAA        mAC        mAG        mAT
  182. //            [C,]    mCA        mCC        mCG        mCT
  183. //            [G,]    mGA        mGC        mGG        mGT
  184. //            [T,]    mTA        mTC        mTG        mTT
  185. //    The general equation for the resulting traces is:
  186. //        R = M O  <==>  R(T,i) =      sum     [ mTS x O(S,i) ]
  187. //                                S in {ACGT}
  188. //    where R is the resulting vector of 4 traces
  189. //          O is the original vector of 4 traces
  190. //          M is the 4x4 ({ACGT}x{ACGT}) matrix whose elements m(i,j) are
  191. //            the cross-term contributions of channel j to channel i
  192. //          S & T are trace identifiers (Source & Target) in {A,C,G,T}
  193. //          i loops over the indices of the trace
  194. //
  195. //    !    The existing traces are discarded and replaced by the transformed data.
  196. //    ========================================|===================================
  197. CTraceFile::error_t
  198. CTraceFile::Transform(xform_mtx& matrix)
  199.     {
  200.     CTrace<trace_t>*    results[NUM_BASES];
  201.  
  202.     // initialize the results array
  203.     for (uint traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  204.         if ( NULL == (results[traceIndex]=new CTrace<trace_t>(numPoints)) )
  205.             // allocate failed... deallocate any already new'd
  206.             {
  207.             do
  208.                 delete results[--traceIndex];
  209.             while(traceIndex);
  210.  
  211.             return memError;
  212.             }
  213.  
  214.     // do the transformations
  215.     uint targetTrace;
  216.     for ( targetTrace = A;targetTrace<NUM_BASES;targetTrace++)
  217.         // targetTrace is the index of the resulting trace
  218.         for (tracepos index=0; index<numPoints; index++)
  219.             // index is the sample point being transformed
  220.             {
  221.             double    sum = 0;
  222.             for (uint sourceTrace=A; sourceTrace<NUM_BASES; sourceTrace++)
  223.                 // sourceTrace is the index of traces contributing to the result
  224.                 sum +=     matrix[targetTrace][sourceTrace]
  225.                                          * (*trace[sourceTrace])[index];
  226.  
  227.             (*results[targetTrace])[index] = (trace_t)sum;
  228.             }
  229.     
  230.     // replace the original traces with the transformed ones.
  231.     for (targetTrace = A; targetTrace<NUM_BASES;targetTrace++)
  232.         {
  233.         delete trace[targetTrace];
  234.         trace[targetTrace] = results[targetTrace];
  235.         trace[targetTrace]->LeftCutoff(leftCutoff);
  236.         trace[targetTrace]->RightCutoff(rightCutoff);
  237.         }
  238.  
  239.     // Calculate the stats for each trace; if mean is <0, then scale by -1
  240.     // orthogonalization will not be affected
  241. //    for (targetTrace = A; targetTrace<NUM_BASES;targetTrace++)
  242. //        {
  243. //        trace[targetTrace]->CalculateStats();
  244. //        if (trace[targetTrace]->Mean()<0)
  245. //            trace[targetTrace]->Scale(-1);
  246. //        }
  247.  
  248.     CalculateStats();
  249.     
  250.     return noError;
  251.     }
  252.  
  253. //    ============================================================================
  254. //    PickPeaks
  255. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  256. //    Calls CTrace<>::PickPeaks for each trace and assimilates the results into a
  257. //    single list in order of occurrence (a merge sort of the peaks).
  258. //    !    Because this is still being designed, error handling is lax, although
  259. //        the skeleton for it exists.
  260. //    ========================================|===================================
  261. CTraceFile::error_t
  262. CTraceFile::PickPeaks(double PeakMeanCoefficient, double ZeroThreshold)
  263.     {
  264.     uint        traceIndex=0;
  265.     //
  266.     // Calculates stats to ensure that each trace has correct min, max, &
  267.     // mean values.
  268.     //
  269.     CalculateStats();
  270.  
  271.     //
  272.     // Generate a list of peaks for each trace
  273.     //
  274.     for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  275.         {
  276.         trace_t MinPeakHeight = (trace_t)(PeakMeanCoefficient
  277.                                              * trace[traceIndex]->Mean());
  278.         if (trace[traceIndex]->PickPeaks(MinPeakHeight,ZeroThreshold)==noError)
  279.             trace[traceIndex]->Peaks().Analyze(*trace[traceIndex]);
  280.         else
  281.             return peakPickError;
  282.         }
  283.  
  284.     return AssimilatePeaks();
  285.     }
  286.  
  287. //    ============================================================================
  288. //    AssimilatePeaks
  289. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  290. //    Assumes that each trace in the tracefile has a list of peaks and collates
  291. //    them into our own list for the tracefile.
  292. //    ========================================|===================================
  293. CTraceFile::error_t
  294. CTraceFile::AssimilatePeaks()
  295.     {
  296.     uint        traceIndex=0;
  297.     //
  298.     // Assimilate the peaks into one CSequence of PeakRec's
  299.     // by popping the 'least' peak from the front of the appropriate
  300.     // list.  When we've done sum total pops, we're done. (This is easier
  301.     // and faster than checking if all lists are empty at every iteration.)
  302.     //
  303.     ulong        sum=0;
  304.     uint        next_least[NUM_BASES]={0};        // next least elem in each trace
  305.  
  306.     for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  307.         sum += trace[traceIndex]->Peaks().Size();
  308.  
  309.     while (sum)
  310.         {
  311.         uint     least_trace = A;
  312.         tracepos least_pos = numPoints;        // just want something that's
  313.                                                 // guaranteed to be larger than
  314.                                                 // our first real hit
  315.  
  316.         // Search for the trace which has the least element (ie. the trace
  317.         // with a peak at an index less than that for peaks of any other trace)
  318.         for (traceIndex=A; traceIndex<NUM_BASES; traceIndex++)
  319.             {
  320.             CPeakList&    ind_peaks = trace[traceIndex]->Peaks();
  321.             if ( (next_least[traceIndex] < ind_peaks.Size()) &&
  322.                  (ind_peaks[next_least[traceIndex]].center < least_pos) )
  323.                 {
  324.                 least_trace = traceIndex;
  325.                 least_pos = ind_peaks[next_least[traceIndex]].center;
  326.                 }
  327.             }
  328.  
  329.         // So now the least peak has been identified.
  330.         // We'll fill in the whichTrace field and append it to the assimilated
  331.         // list of peaks
  332.         PeakRec pr = trace[least_trace]->Peaks()[next_least[least_trace]];
  333.         pr.whichTrace = (base_t)least_trace;
  334.         peaks.Append(pr);
  335.  
  336.         // Set the next_least index to the successor of this victim and
  337.         // decrement the sum counter... 1 down and sum [sic] more to go.  :-|
  338.         next_least[least_trace]++;
  339.         sum--;
  340.         }
  341.     return noError;
  342.     }
  343.  
  344. //    ============================================================================
  345. //    PrunePeaks
  346. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  347. //    Simply step through the sequence of peaks and in the case where adjacent
  348. //    peaks are less than minSeparation apart, remove the peak with the lesser
  349. //    probability.
  350. //    A list of the pairwise peak contentions and their resolution will be output
  351. //    if an output stream is provided.
  352. //
  353. //    ALGORITHM:
  354. //    Given a set of peak records indexed 0..s  with fields:
  355. //        center - position in trace of the center of the peak
  356. //        PH - probability of the peak
  357. //    and a minimum separation, minSep, prune the list by removing the less
  358. //    probable peak when adjacent peak centers are within minSep.
  359. //    let peakA be the index of the first peak in the comparison
  360. //    let peakA=0
  361. //    while peakA<=s-2 do
  362. //        let peakB = peakA + 1
  363. //        if peakB.center - peakA.center <= minSep then
  364. //            if peakA.PH < peakB.PH then
  365. //                remove peakA
  366. //            else
  367. //                remove peakB
  368. //            fi
  369. //        else
  370. //            increment peakA
  371. //        fi
  372. //    elihw
  373. //
  374. //    !    This is intended to be only a first cut at a pruning algorithm.  It has
  375. //        many problems, including the possibility that a series of peaks all
  376. //        separated by less than minSeparation between adjacent members will
  377. //        result in only 1 of the members being selected, when perhaps more than
  378. //        one peak should result from the collection.
  379. //    ========================================|===================================
  380. CTraceFile::error_t
  381. CTraceFile::PrunePeaks(
  382.     tracepos    minSeparation,
  383.     ostream*    os)
  384.     {
  385.     ulong peakNo=0;
  386.  
  387.     while(peakNo<=peaks.Size()-2)
  388.         if (peaks[peakNo+1].center-peaks[peakNo].center < minSeparation)
  389.             {
  390.             // victim is the peak which is less probable [ie. lesser P(H|D)]
  391.             ulong victim = (peaks[peakNo+1].PHD < peaks[peakNo].PHD ?
  392.                             peakNo+1:
  393.                             peakNo);
  394.             base_t victimTrace = peaks[victim].whichTrace;
  395.  
  396.             if (os)
  397.                 {
  398.                 *os << setfill('-') << setw(80) << "" << setfill(' ') << endl;
  399.                 *os << Tag(victim==peakNo)   << peaks[peakNo];
  400.                 *os << Tag(victim==peakNo+1) << peaks[peakNo+1];
  401.                 }
  402.  
  403.             // remove the peak record from the trace's peak list
  404.             (*trace[victimTrace]).Peaks().RemovePeak(peaks[victim].center);
  405.  
  406.             // and remove the peak from the tracefile's peak list
  407.             peaks.Remove(victim);
  408.             }
  409.         else
  410.             peakNo++;
  411.         // end while
  412.  
  413.     if (os)
  414.         *os << setfill('-') << setw(80) << "" << setfill(' ') << endl;
  415.  
  416.     return noError;
  417.     }
  418.  
  419. ostream&
  420. operator<<(ostream& os, CTraceFile& ctf)
  421.     {
  422.     tracepos length=ctf.rightCutoff - ctf.leftCutoff + 1;
  423.           os    << "filename:\t"     << ctf.filename                << endl
  424.         << "native format:\t" << FileFormatAbbr(ctf.nativeFormat)
  425.             << " [" << FileFormatDesc(ctf.nativeFormat) << "]" << endl
  426.         << "strand:\t\t"    << (short)ctf.whichStrand            << endl
  427.         << "total size:\t"    << ctf.numPoints            << endl
  428.         << "window:\t\t"
  429.             << "[" << ctf.leftCutoff << "," << ctf.rightCutoff << "]"
  430.             << "  (" << length << " point" << Plural(length) << ")" << endl
  431.         << "min/max/mean:\t" << ctf.min << "/"
  432.                             << ctf.max << "/"
  433.                             << ctf.mean << "  (within window)" << endl
  434.         << "sequence len:\t"<< ctf.numBases            << endl
  435.         << "comments" << endl << "--------" << endl << ctf.comments << endl
  436.         << "peaks picked:\t"<< ctf.peaks.Size()            << endl;
  437.  
  438.     if (ctf.peaks.Size() != 0)
  439.         {
  440.         os << setw(3) << "#" << "  " << PeakRecHeader;
  441.         for(uint i = 0; i<ctf.peaks.Size(); i++)
  442.             os << setw(3) << i+1 << "  " << ctf.peaks[i];
  443.         }
  444.  
  445.     return os;
  446.     }
  447.  
  448. //    ============================================================================
  449. //    ###        # #######
  450. //     #        #  #     #
  451. //     #       #   #     #
  452. //     #      #    #     #
  453. //     #     #     #     #
  454. //     #    #      #     #
  455. //    ###  #       #######
  456. //    
  457. //    General I/O Methods
  458. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  459. //    public:
  460. //    Read(const char*, format) - Reads from named file in specified format, or
  461. //        inferred format if not specified.
  462. //    Write(const char*, format) - Writes to named file in specified format, or
  463. //        native format if not specified.
  464. //
  465. //    private:
  466. //    Read(FILE*, format) - Reads from FILE* (opened with at least "rb"
  467. //        permissions in specified format (mandatory).  This just dispatches to
  468. //        the appropriate file format reader (ie. ReadSCF).
  469. //    Write(FILE*, format) - Writes to FILE* (opened with at least "wb"
  470. //        permissions) in specified format (mandatory).  This just dispatches to
  471. //        the appropriate file format writer (ie. WriteSCF).
  472. //
  473. //    NOTES
  474. //    *    File I/O must be handled through the const char* reader and writer.
  475. //        All other methods are declared private.
  476. //    *    The private methods do not set the class's error flag; this allows the
  477. //        calling function (ie. the public methods) to handle the error as
  478. //        appropriate, or to set the error flag if desired.  (ie. we may wish to
  479. //        retry with a different format.
  480. //    ?    Currently, if an error occurs during reading, it's possible to return
  481. //        (with error) an incompletely-read trace.  I think that this is not
  482. //        desired, but I'm still weighing the merits.
  483. //    ========================================|===================================
  484.  
  485. //    ============================================================================
  486. //    Read(const char*, format_t)
  487. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  488. //    Takes a filename and optional format.  If the format is not specified, I
  489. //    attempt to determine it by FileFormat (see FileFormat.H).  If it still can't
  490. //    be determined, we're sunk; if it can be determined, we attempt to read it.
  491. //
  492. //    NOTES
  493. //    ?    should the trace be complete deallocated if it fails?
  494. //    ========================================|===================================
  495. CTraceFile::error_t
  496. CTraceFile::Read(const char* fn, format_t fmt)
  497.     {
  498.     error_t        err;
  499.     FILE*        fp;
  500.  
  501.     if (unknown == fmt)                        // user didn't tell us format... try
  502.         {                                    //   to determine it for ourselves
  503.         fmt = FileFormat(fn);
  504.         if (unknown == fmt)                    // drat! still unknown
  505.             return (error = unkFmtError);
  506.         }
  507.  
  508.     if (fp = fopen(fn,"rb"))
  509.         {
  510.         // file opened successfully
  511.         err = Read(fp,fmt);
  512.         fclose(fp);
  513.         }
  514.     else
  515.         {
  516.         err = fileDoesntExistError;
  517.         }
  518.  
  519.     if (!err)
  520.         {
  521.         filename = new char[strlen(fn)];
  522.         strcpy(filename,fn);
  523.         nativeFormat = fmt;
  524.  
  525.         CalculateStats();
  526.         }
  527.     
  528.     return (error = err);
  529.     }
  530.  
  531. //    ============================================================================
  532. //    Write(const char*, format_t)
  533. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  534. //    Takes a filename and optional format.  If the format is not specified, the
  535. //    native format is assumed.  If we're writing in format other than the one
  536. //    used to read the tracefile, the comment field is duplicated and a temporary
  537. //    version is constructed which annotates the conversion process.
  538. //
  539. //    NOTES
  540. //    !    I don't check to see if the file exists.  I should do this...
  541. //    ========================================|===================================
  542. CTraceFile::error_t
  543. CTraceFile::Write(const char* fn, format_t fmt)
  544.     {
  545.     error_t        err;
  546.     FILE*        fp;
  547.     char*        tempComment=NULL;            // temporary comment
  548.     char*        origComment=comments;        // saved pointer to original comment
  549.  
  550.     if (DEBUG)    cerr << "Write(char*,format_t): Entering" << endl;
  551.  
  552.     if ( 0 == numPoints )
  553.         return (error = emptyError);
  554.     
  555.     if (fmt == unknown)                        // user didn't tell us format...
  556.         fmt = nativeFormat;                    //   assume nativeFormat
  557.     else
  558.         if (fmt != nativeFormat)
  559.             // user specified a format different from native... make a comment
  560.             // that this was translated by CTraceFile.
  561.             {
  562.             char buf[1000] = {0};
  563.  
  564.             strcat(buf,"conv = by CTraceFile ");
  565.             strcat(buf,version);
  566.             strcat(buf," from ");
  567.             strcat(buf,FileFormatDesc(nativeFormat));
  568.             strcat(buf," (");
  569.             strcat(buf,FileFormatAbbr(nativeFormat));
  570.             strcat(buf,")");
  571.             strcat(buf," to ");
  572.             strcat(buf,FileFormatDesc(fmt));
  573.             strcat(buf," (");
  574.             strcat(buf,FileFormatAbbr(fmt));
  575.             strcat(buf,")");
  576.             strcat(buf," on ");
  577.             time_t tm = time(NULL);
  578. //            char* tmsp = ctime(&tm);
  579. //            strcat(buf,tmsp);
  580.             struct tm* tp=localtime(&tm);
  581.             strftime(buf+strlen(buf),1000-strlen(buf),"%Y-%b-%d %H:%M%Z\n",tp);
  582.  
  583. //            This is the original conversion comment; I've opted for a style
  584. //            that more closely matches makeSCF's for consistency.
  585. //            strcat(buf,tmsp);
  586. //            strcpy(buf+strlen(buf)-1,"- tracefile converted from ");
  587. //            strcat(buf,FileFormatDesc(nativeFormat));
  588. //            strcat(buf," to ");
  589. //            strcat(buf,FileFormatDesc(fmt));
  590. //            strcat(buf," by CTraceFile.\n");
  591.                         
  592.             // make space for the temporary comment
  593.             tempComment = new char[strlen(comments)+strlen(buf)+1];
  594.             if (tempComment == NULL)
  595.                 return memError;
  596.             strcpy(tempComment,comments);    // copy old comment to new comment
  597.             strcat(tempComment,buf);        // add the conversion comment
  598.  
  599.             comments = tempComment;            // we'll write, then swap them back
  600.             }
  601.  
  602.     if (DEBUG)    cerr << "Write(char*,format_t): Writing" << endl;
  603.  
  604.     if ((fp = fopen(fn,"wb")) != NULL)
  605.         {
  606.         // file opened successfully
  607.         if ((err = Write(fp,fmt)) != noError)
  608.             {
  609.             fclose(fp);
  610.             remove(fn);                        // half-complete files are useless
  611.             }
  612.         else
  613.             fclose(fp);                        // write succeeded
  614.         }
  615.     else
  616.         err = ioError;
  617.  
  618.     if (DEBUG)    cerr << "Write(char*,format_t): Exiting" << endl;
  619.  
  620.     if (tempComment != NULL)                // this write involved a conversion
  621.         {                                    // we need to restore the original
  622.         comments=origComment;                // comment and kill the space for
  623.         delete tempComment;                    // the temporary
  624.         }
  625.  
  626.     return (error = err);
  627.     }
  628.  
  629. //    ============================================================================
  630. //    Read(FILE*, format_t) [private]
  631. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  632. //    This is the dispatcher to the appropriate format reader.  It is called by
  633. //    Read(const char*, format_t), but not called directly.
  634. //    ========================================|===================================
  635. CTraceFile::error_t
  636. CTraceFile::Read(FILE* fp, format_t fmt)
  637.     {
  638.     if (DEBUG)    cerr << "Read(FILE*,format_t): Entering; fmt = "<< (short)fmt<< endl;
  639.  
  640.     if (NULL == fp)
  641.         return ioError; 
  642.  
  643.     // else...
  644.     error_t    err = noError;
  645.     switch (fmt)
  646.         {
  647.         case unknown:
  648.             err = unkFmtError;
  649.             break;
  650.         case SCF:
  651.             err = ReadSCF(fp);
  652.             break;
  653.         case ABI0:
  654.             err = ReadABI(fp,0);
  655.             break;
  656.         case ABI1:
  657.             err = ReadABI(fp,1);
  658.             break;
  659.         case ABI:
  660.             err = ReadABI(fp,2);
  661.             break;
  662.         case ALF:
  663.             err = ReadALF(fp);
  664.             break;
  665.         default:
  666.             err = fmtNotSuppError;
  667.         }
  668.  
  669.     if (DEBUG)    cerr << "Read(FILE*,format_t): Exiting" << endl;
  670.     return err;
  671.     }
  672.  
  673. //    ============================================================================
  674. //    Write(FILE*, format_t) [private]
  675. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  676. //    This is the dispatcher to the appropriate format writer.  It will most
  677. //    likely be called by Write(const char*, format_t), but not called directly.
  678. //    ========================================|===================================
  679. CTraceFile::error_t
  680. CTraceFile::Write(FILE* fp, format_t fmt)
  681.     {
  682.     error_t    err = noError;
  683.  
  684.     if (NULL == fp)
  685.         return ioError;
  686.  
  687.     // else...
  688.     if (DEBUG)    cerr << "Write(FILE*,format_t): Entering" << endl;
  689.  
  690.     switch (fmt)
  691.         {
  692.         case unknown:
  693.             err = unkFmtError;
  694.             break;
  695.         case SCF:
  696.             err = WriteSCF(fp);
  697.             break;
  698.         case ABI:    
  699.             err = WriteABI(fp);
  700.             break;
  701.         case ALF:
  702.             err = WriteALF(fp);
  703.             break;
  704.         default:
  705.             err = fmtNotSuppError;
  706.         }
  707.  
  708.     if (DEBUG)    cerr << "Write(FILE*,format_t): Exiting" << endl;
  709.     return err;
  710.     }
  711.  
  712.  
  713. //    ============================================================================
  714. //     #####   #####  #######           ###         # #######
  715. //    #     # #     # #                  #         #  #     #
  716. //    #       #       #                  #        #   #     #
  717. //     #####  #       #####              #       #    #     #
  718. //          # #       #                  #      #     #     #
  719. //    #     # #     # #                  #     #      #     #
  720. //     #####   #####  #                 ###   #       #######
  721. //    
  722. //    SCF (Standard Chromatogram Format) I/O Methods
  723. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  724. //    File format inferred from ted source (see CTraceFile.H).
  725. //
  726. //    SCF_File = 
  727. //        offset            description
  728. //        ---------------    --------------------------------
  729. //        0                <SCF_Header>
  730. //        trace_offset    <SCF_Sample> x num_samples
  731. //        sequence_offset    <SCF_Bases> x num_bases
  732. //        comment_offset    comment of length comment_length
  733. //    ========================================|===================================
  734. const long SCF_MAGIC_COOKIE = ((('.'<<8)+'s'<<8)+'c'<<8)+'f'; 
  735.  
  736. struct    SCF_Header
  737.     {
  738.     long    magic_cookie;                    // = SCF_MAGIC_COOKIE
  739.     long    num_points;                        // number of trace points (per base)
  740.     long    trace_offset;                    // position of trace data in file
  741.     long    num_bases;                        // number of bases in sequence
  742.     long    left_cutoff;                    // unreadable < left cutoff
  743.     long    right_cutoff;                    // unreadable > right cutoff
  744.     long    sequence_offset;                // position of sequence data in file
  745.     long    comment_length;                    // strlen of comment
  746.     long    comment_offset;                    // position of comment in file
  747.     long    spare[23];
  748.  
  749.     friend size_t fread(FILE* fp, SCF_Header& header)  
  750.         { return ::fread(&header,sizeof(SCF_Header),1,fp); }
  751.     friend size_t fwrite(FILE* fp, const SCF_Header& header)  
  752.         { return ::fwrite(&header,sizeof(SCF_Header),1,fp); }
  753.     };
  754.  
  755. struct    SCF_Sample                            // note that the values for each
  756.     {                                        // trace are interleaved. ie
  757.     byte    sample_A;                        // A1-C1-G1-T1-A2-C2-G2-...-Gn-Tn
  758.     byte    sample_C;                        // where Ni = value of point i of
  759.     byte    sample_G;                        // trace N.
  760.     byte    sample_T;
  761.  
  762.     friend size_t    fread(FILE* fp, SCF_Sample& sample)  
  763.         { return ::fread(&sample,sizeof(SCF_Sample),1,fp); }
  764.     friend size_t    fwrite(FILE* fp, const SCF_Sample& sample)
  765.         { return ::fwrite(&sample,sizeof(SCF_Sample),1,fp); }
  766.     };
  767.  
  768. struct    SCF_Base
  769.     {
  770.     long    trace_index;                    // location of base as called
  771.     byte    prob_A;                            // probability that base = A
  772.     byte    prob_C;                            //                           C
  773.     byte    prob_G;                            //                           G
  774.     byte    prob_T;                            //                           T
  775.     char    base;                            // base that was chosen
  776.     byte    _unused1;                        // unused?
  777.     byte    _unused2;                        // unused?
  778.     byte    _unused3;                        // unused?
  779.  
  780.     friend size_t fread(FILE* fp, SCF_Base& base)  
  781.         { return ::fread(&base,sizeof(SCF_Base),1,fp); }
  782.     friend size_t fwrite(FILE* fp, const SCF_Base& base) 
  783.         { return ::fwrite(&base,sizeof(SCF_Base),1,fp); }
  784.     };
  785.  
  786. //    ============================================================================
  787. //    ReadSCF [private]
  788. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  789. //    Read a SCF file; called exclusively by Read(FILE*,format_t).
  790. //    ========================================|===================================
  791. CTraceFile::error_t
  792. CTraceFile::ReadSCF(FILE* fp)
  793.     {
  794.     uint        position;
  795.  
  796.     if (ferror(fp))
  797.         return ioError;
  798.  
  799.     //
  800.     // Read header
  801.     //
  802.     SCF_Header    header;
  803.  
  804.     rewind(fp);                                // take it from the top...
  805.     fread(fp,header);                        // Read the header
  806.     if (ferror(fp))                            // check for error
  807.         return ioError;
  808.  
  809.     nativeFormat = SCF;                        // set our instance's data members
  810.     whichStrand = top;
  811.  
  812.  
  813.     //
  814.     // Allocate space for traces, sequence, etc.
  815.     //   this is all-or-none: we either want all our required space or none at
  816.     //   all.
  817.     error_t err = Allocate((size_t)header.num_points,\
  818.                            (size_t)header.num_bases,\
  819.                            (size_t)header.comment_length);
  820.     if (noError != err)
  821.         return err;
  822.  
  823.     //
  824.     // read the traces
  825.     //
  826.     fseek(fp,header.trace_offset,SEEK_SET);
  827.     for (position=0; position<numPoints; position++)
  828.         {
  829.         SCF_Sample    sample;
  830.         fread(fp,sample);
  831.         (*trace[A])[position]=sample.sample_A;
  832.         (*trace[C])[position]=sample.sample_C;
  833.         (*trace[G])[position]=sample.sample_G;
  834.         (*trace[T])[position]=sample.sample_T;
  835.         }
  836.     if (ferror(fp))
  837.         return ioError;
  838.  
  839.     //
  840.     // read the sequence & base positions
  841.     //
  842.     fseek(fp,header.sequence_offset,SEEK_SET);
  843.     for (uint traceIndex=0; traceIndex<numBases; traceIndex++)
  844.         {
  845.         SCF_Base    base;
  846.         fread(fp,base);
  847.         sequence[traceIndex]=base.base;
  848.         basePositions[traceIndex]=(tracepos)base.trace_index;
  849.         }
  850.     if (ferror(fp))
  851.         return ioError;
  852.  
  853.     //
  854.     // read the comment
  855.     //
  856.     fseek(fp,header.comment_offset,SEEK_SET);
  857.     fread(comments,(size_t)header.comment_length,1,fp);
  858.     if (ferror(fp))
  859.         return ioError;
  860.  
  861.     return noError;
  862.     }
  863.  
  864. //    ============================================================================
  865. //    WriteSCF [private]
  866. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  867. //    Write a SCF file; called exclusively by Write(FILE*,format_t).
  868. //    ========================================|===================================
  869. CTraceFile::error_t
  870. CTraceFile::WriteSCF(FILE* fp)
  871.     {
  872.     SCF_Header    header;
  873.     uint        position;
  874.  
  875.     if (DEBUG)    cerr << "WriteSCF: Entering" << endl;
  876.  
  877.     //
  878.     // prepare the header and calculate file offsets for trace information
  879.     //
  880.     header.magic_cookie        = SCF_MAGIC_COOKIE;
  881.     header.num_points        = numPoints;
  882.     header.num_bases        = numBases;
  883.     header.left_cutoff        = leftCutoff;
  884.     header.right_cutoff        = rightCutoff;
  885.     header.comment_length    = strlen(comments);
  886.     header.trace_offset        = sizeof(SCF_Header);
  887.     header.sequence_offset    = header.trace_offset+sizeof(SCF_Sample)*numPoints;
  888.     header.comment_offset    = header.sequence_offset+sizeof(SCF_Base)*numBases;
  889.  
  890.     if (DEBUG)    cerr << "WriteSCF: header" << endl;
  891.  
  892.     //
  893.     // write the header
  894.     //
  895.     rewind(fp);                                // take it from the top...
  896.     fwrite(fp, header);                        // write the header
  897.     if (ferror(fp))                            // check for error
  898.         return ioError;
  899.  
  900.     if (DEBUG)    cerr << "WriteSCF: traces" << endl;
  901.  
  902.     //
  903.     // write the traces
  904.     //
  905.     fseek(fp,header.trace_offset,SEEK_SET);    // move to traces
  906.     for (position=0; position<numPoints; position++)
  907.         {
  908.         SCF_Sample    sample;
  909.         sample.sample_A=(byte)(*trace[A])[position];
  910.         sample.sample_C=(byte)(*trace[C])[position];
  911.         sample.sample_G=(byte)(*trace[G])[position];
  912.         sample.sample_T=(byte)(*trace[T])[position];
  913.         fwrite(fp,sample);
  914.         }
  915.     if (ferror(fp))
  916.         return ioError;
  917.  
  918.     if (DEBUG)    cerr << "WriteSCF: sequence" << endl;
  919.  
  920.     //
  921.     // write the sequence & base positions
  922.     //
  923.     fseek(fp,header.sequence_offset,SEEK_SET);    // move to traces
  924.     for (position=0; position<numBases; position++)
  925.         {
  926.         SCF_Base    base;
  927.         base.base=sequence[position];
  928.         base.trace_index=basePositions[position];
  929.         fwrite(fp,base);
  930.         }
  931.     if (ferror(fp))
  932.         return ioError;
  933.  
  934.     if (DEBUG)    cerr << "WriteSCF: comment" << endl;
  935.  
  936.     //
  937.     // write the comment
  938.     //
  939.     fseek(fp,header.comment_offset,SEEK_SET);
  940.     fwrite(comments,(size_t)header.comment_length,1,fp);
  941.     if (ferror(fp))
  942.         return ioError;
  943.  
  944.     if (DEBUG)    cerr << "WriteSCF: Exiting" << endl;
  945.  
  946.     return noError;
  947.     }
  948.  
  949. //    ============================================================================
  950. //       #    ######    ###             ###         # #######
  951. //      # #   #     #    #               #         #  #     #
  952. //     #   #  #     #    #               #        #   #     #
  953. //    #     # ######     #               #       #    #     #
  954. //    ####### #     #    #               #      #     #     #
  955. //    #     # #     #    #               #     #      #     #
  956. //    #     # ######    ###             ###   #       #######
  957. //
  958. //    ABI format I/O methods
  959. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  960. //    These routines are copied almost verbatim from seqIOABI.[hc] written by
  961. //    LaDeana Hillier, although I've made extensive modifications to certain
  962. //    sections.
  963. //    ========================================|===================================
  964. typedef ushort int2;                        // Two byte short assumption
  965. typedef ulong  int4;                        // Four byte long assumption
  966.  
  967. const long DATAtag            = ((long) ((((('D'<<8)+'A')<<8)+'T')<<8)+'A');
  968. const long BASEtag            = ((long) ((((('P'<<8)+'B')<<8)+'A')<<8)+'S');
  969. const long BASEPOStag        = ((long) ((((('P'<<8)+'L')<<8)+'O')<<8)+'C');
  970. const long SPACINGtag        = ((long) ((((('S'<<8)+'P')<<8)+'A')<<8)+'C');
  971. const long SIGNALtag        = ((long) ((((('S'<<8)+'/')<<8)+'N')<<8)+'%');
  972. const long FWO_tag            = ((long) ((((('F'<<8)+'W')<<8)+'O')<<8)+'_');
  973. const long MCHNtag            = ((long) ((((('M'<<8)+'C')<<8)+'H')<<8)+'N');
  974. const long PDMFtag            = ((long) ((((('P'<<8)+'D')<<8)+'M')<<8)+'F');
  975. const long SMPLtag            = ((long) ((((('S'<<8)+'M')<<8)+'P')<<8)+'L');
  976. const long PPOStag            = ((long) ((((('P'<<8)+'P')<<8)+'O')<<8)+'S');
  977. const long PRIMERtag        = ((long) ((((('P'<<8)+'P')<<8)+'O')<<8)+'S');
  978.  
  979. const size_t INDEXPO        = 26;            // offset of offset to index
  980. const int    INDEX_ENTRY_LENGTH= 28;
  981.  
  982. inline
  983. short
  984. baseIndex(char B)
  985.     {
  986.     return ((B)=='C'?0:(B)=='A'?1:(B)=='G'?2:3);
  987.     }
  988.  
  989. //    ============================================================================
  990. //    readABIInt2
  991. //    ========================================|===================================
  992. static
  993. bool
  994. readABIInt2(
  995.     FILE    *fp,
  996.     int2    *i2)
  997.     {
  998.     unsigned char buf[sizeof(int2)];
  999.  
  1000.     if (fread(buf, sizeof(buf), 1, fp) != 1) return (FALSE);
  1001.     *i2 = (int2)
  1002.         (((unsigned short)buf[1]) +
  1003.          ((unsigned short)buf[0]<<8));
  1004.     return (TRUE);
  1005.     }
  1006.  
  1007. //    ============================================================================
  1008. //    readABIInt4
  1009. //    ========================================|===================================
  1010. static
  1011. bool
  1012. readABIInt4(
  1013.     FILE    *fp,
  1014.     int4    *i4)
  1015.     {
  1016.     unsigned char buf[sizeof(int4)];
  1017.  
  1018.     if (fread(buf, sizeof(buf), 1, fp) != 1) return (FALSE);
  1019.     *i4 = (int4)
  1020.         (((unsigned long)buf[3]) +
  1021.          ((unsigned long)buf[2]<<8) +
  1022.          ((unsigned long)buf[1]<<16) +
  1023.          ((unsigned long)buf[0]<<24));
  1024.     return (TRUE);
  1025.     }
  1026.  
  1027. //    ============================================================================
  1028. //    getIndexEntryLW
  1029. //    -  -  -     -    -  -  -     -    -  -  -     -    -  -  -     -    -  -  -     -    -  -  -     -    -  -
  1030. //    From the ABI results file connected to `fp' whose index starts at byte
  1031. //    offset `indexO', return in `val' the `lw'th long word from the `count'th
  1032. //    entry labeled `label'. The result indicates success.
  1033. //    ========================================|===================================
  1034. bool
  1035. getIndexEntryLW (
  1036.     FILE*    fp,
  1037.     long    indexO,
  1038.     long    label,
  1039.     long    count,
  1040.     int        lw,
  1041.     int4*    val)
  1042.     {
  1043.     int        entryNum=-1;
  1044.     int        i;
  1045.     int4    entryLabel, entryLw1;
  1046.  
  1047. //    if (DEBUG)
  1048. //cerr<<"getIndexEntryLW(label="<<label<<",count="<<count<<",lw="<<lw<<")"<< endl;
  1049.  
  1050.     do
  1051.         {
  1052.         entryNum++;
  1053.         if (fseek(fp, indexO+(entryNum*INDEX_ENTRY_LENGTH), SEEK_SET) != 0)
  1054.             return FALSE;
  1055.         if (!readABIInt4(fp, &entryLabel))
  1056.             return FALSE;
  1057.         if (!readABIInt4(fp, &entryLw1))
  1058.             return FALSE;
  1059.         } while (!(entryLabel == label && entryLw1 == count));
  1060.  
  1061.     for (i=2; i<=lw; i++)
  1062.         if (!readABIInt4(fp, val)) return FALSE;
  1063.  
  1064.     return TRUE;
  1065.     }
  1066.  
  1067. //    ============================================================================
  1068. //    getIndexEntryW
  1069. //    -  -  -     -    -  -  -     -    -  -  -     -    -  -  -     -    -  -  -     -    -  -  -     -    -  -
  1070. //    From the ALF results file connected to `fp' whose index starts at byte
  1071. //    offset `indexO', return in `val' the `lw'th     word (int2) from the entry
  1072. //    labeled `label'. The result indicates success.
  1073. //    ========================================|===================================
  1074. bool
  1075. getIndexEntryW (
  1076.     FILE*    fp,
  1077.     long    indexO,
  1078.     long    label,
  1079.     int        lw,
  1080.     int2*    val)
  1081.     {
  1082.     int        entryNum=-1;
  1083.     int        i;
  1084.     int4    entryLabel;
  1085.     int4    jval;
  1086.  
  1087.     do
  1088.         {
  1089.         entryNum++;
  1090.         if (fseek(fp, indexO+(entryNum*INDEX_ENTRY_LENGTH), SEEK_SET) != 0)
  1091.             return FALSE;
  1092.         if (!readABIInt4(fp, &entryLabel))
  1093.             return FALSE;
  1094.         }
  1095.         while (entryLabel != label);
  1096.  
  1097.  
  1098.     for (i=2; i<lw; i++)
  1099.         if (!readABIInt4(fp, &jval)) return FALSE;
  1100.     
  1101.     if (!readABIInt2(fp, val)) return FALSE;
  1102.  
  1103.     return TRUE;
  1104.     }
  1105.  
  1106. //    ============================================================================
  1107. //    ReadABI [private]
  1108. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  1109. //    Read a ABI file; called exclusively by Read(FILE*,format_t).
  1110. //    Read the ABI format sequence with name `fn' into `seq'.
  1111. //
  1112. //    This code is a rehash of seqIOABI.c, significantly modified for use here.
  1113. //    bottom strand designations are completely ignored here.
  1114. //    ========================================|===================================
  1115. CTraceFile::error_t
  1116. CTraceFile::ReadABI(FILE* fp, short whichSet)
  1117.     {
  1118.     if (DEBUG)    cerr << "ReadABI: Entering; set = " << whichSet << endl;
  1119.  
  1120.     const size_t COMMENT_LEN    = 2048;        // max comment length;
  1121.     const short TRACE_INDEX        = NUM_BASES*whichSet+1;
  1122.                                     // index to offset of 1st trace in set
  1123.  
  1124.     error_t    err = noError;
  1125.     char    line[128];
  1126.     int4    ppos;
  1127.  
  1128.     bool    bottom = FALSE;
  1129.     char*    enzyme = NULL;
  1130.  
  1131.     uint    whichTrace;
  1132.     ulong    numPoints;
  1133.     ulong    numBases;
  1134. //    int2    primerPos;
  1135.     int2    traceValue;
  1136.  
  1137.     int4    fwo_;                            // base -> lane mapping
  1138.     int4    signalO;
  1139.     int4    indexO;                            // index
  1140.     int4    baseO;                            // bases
  1141.     int4    basePosO;                        // base positions
  1142.     int4    dataCO;                            // C trace
  1143.     int4    dataAO;                            // A trace
  1144.     int4    dataGO;                            // G trace
  1145.     int4    dataTO;                            // T trace
  1146.     int4    MCHN_O;                            // machine name
  1147.     int4    PDMF_O;                            // dye primer guff
  1148.     int4    SMPL_O;                            // sample name
  1149.  
  1150.     int        i;
  1151.  
  1152.     if (DEBUG)    cerr << "ReadABI: 1.0" << endl;
  1153.     // get the index offset
  1154.     if ((fseek(fp, INDEXPO, SEEK_SET) != 0) || (!readABIInt4(fp, &indexO)))
  1155.         return ioError;
  1156.  
  1157.     if (DEBUG)    cerr << "ReadABI: 1.1" << endl;
  1158.     // get the number of points
  1159.     if (!getIndexEntryLW(fp,indexO,DATAtag,TRACE_INDEX,3,&numPoints))
  1160.         return ioError;
  1161.     
  1162.     // get the number of bases
  1163.     if (DEBUG)    cerr << "ReadABI: 1.2" << endl;
  1164.     if (!getIndexEntryLW(fp,indexO,BASEtag,1,3,&numBases))
  1165.         return ioError;
  1166.  
  1167.     //    not needed; in fact it's apparently superfluous and incorrect when
  1168.     //    compared with similar code to do the same thing, below.  -- Reece
  1169.     #if (FALSE)
  1170.     // get Primer Position...where the ABI file found the start of
  1171.     if (DEBUG)    cerr << "ReadABI: 1.3" << endl;
  1172.     if (!getIndexEntryW(fp,indexO,PPOStag,6,&primerPos))
  1173.         return ioError;
  1174.     #endif
  1175.  
  1176.     if (DEBUG)    cerr << "ReadABI: 1.4" << endl;
  1177.     // get bases offset
  1178.     if (!getIndexEntryLW(fp,indexO,BASEtag,1,5,&baseO))
  1179.         return ioError;
  1180.  
  1181.     if (DEBUG)    cerr << "ReadABI: 1.5" << endl;
  1182.     // get base positions offset
  1183.     if (!getIndexEntryLW(fp,indexO,BASEPOStag,1,5,&basePosO))
  1184.         return ioError;
  1185.  
  1186.     if (DEBUG)    cerr << "ReadABI: 1.6" << endl;
  1187.     // The order of the DATA fields is determined by the field FWO_
  1188.     // Juggle around with data pointers to get it right
  1189.     int4*    dataxO[NUM_BASES];
  1190.     dataxO[0] = &dataCO;
  1191.     dataxO[1] = &dataAO;
  1192.     dataxO[2] = &dataGO;
  1193.     dataxO[3] = &dataTO;
  1194.  
  1195.     // Get the Freak World Out (FWO?) field ...
  1196.     if (!getIndexEntryLW(fp,indexO,FWO_tag,1,5,&fwo_))
  1197.         return ioError;
  1198.     //Get the positions of the four traces
  1199.     if (!(
  1200.         getIndexEntryLW(fp,indexO,DATAtag,TRACE_INDEX,
  1201.                         5,dataxO[baseIndex((char)(fwo_>>24&BYTE[0]))]) &&
  1202.         getIndexEntryLW(fp,indexO,DATAtag,TRACE_INDEX+1,
  1203.                         5,dataxO[baseIndex((char)(fwo_>>16&BYTE[0]))]) &&
  1204.         getIndexEntryLW(fp,indexO,DATAtag,TRACE_INDEX+2,
  1205.                         5,dataxO[baseIndex((char)(fwo_>>8&BYTE[0]))]) &&
  1206.         getIndexEntryLW(fp,indexO,DATAtag,TRACE_INDEX+3,
  1207.                         5,dataxO[baseIndex((char)(fwo_&BYTE[0]))]) ))
  1208.         return ioError;
  1209.  
  1210.  
  1211.     //
  1212.     // Offsets for critical data have been obtained.  Offsets for non-critical
  1213.     // info will be determined later, mostly for comment field.
  1214.     // Now we're reade to get the critical information.
  1215.     //
  1216.     if (DEBUG)    cerr << "ReadABI: 2-allocate" << endl;
  1217.  
  1218.     //
  1219.     // Allocate space for data, bases, and comments
  1220.     //
  1221.     if (err = Allocate((size_t)numPoints,(size_t)numBases,COMMENT_LEN))
  1222.         return err;
  1223.  
  1224.     //
  1225.     // Read in the traces using the offsets above
  1226.     //
  1227.     int4    dataOffset[NUM_BASES];
  1228.             dataOffset[A]=dataAO;    dataOffset[C]=dataCO;
  1229.             dataOffset[G]=dataGO;    dataOffset[T]=dataTO;
  1230.  
  1231.     for (whichTrace=A;whichTrace<=T;whichTrace++)
  1232.         {
  1233.         if (DEBUG)    cerr << "ReadABI: 2" << DNA_BASES[whichTrace]
  1234.                         << "; offset = " << dataOffset[whichTrace] << endl;
  1235.  
  1236.         if (fseek(fp, dataOffset[whichTrace], SEEK_SET))
  1237.             return ioError;
  1238.         for (i=0;i<numPoints;i++)
  1239.             if (readABIInt2(fp, &traceValue))
  1240.                 (*trace[whichTrace])[i] = traceValue;
  1241.             else
  1242.                 return ioError;
  1243.         }
  1244.  
  1245.     //
  1246.     // read in the sequence
  1247.     //
  1248.     if (DEBUG)    cerr << "ReadABI: 3" << endl;
  1249.     if (fseek(fp, baseO, SEEK_SET))
  1250.         return ioError;
  1251.     for (i=0;i<numBases;i++)
  1252.         {
  1253.         int        ch;
  1254.         if ((ch = fgetc(fp)) == EOF)
  1255.             return ioError;
  1256.         sequence[i] = (ch == 'N') ? '-' : ch;
  1257.         }
  1258.  
  1259.     //
  1260.     // read base positions
  1261.     //
  1262.     if (DEBUG)    cerr << "ReadABI: 4" << endl;
  1263.     for (i=0;i<numBases;i++)
  1264.         if (!readABIInt2(fp, (int2*)&basePositions[i]))
  1265.             return ioError;
  1266.  
  1267.  
  1268.     //
  1269.     // Gather useful information
  1270.     // such as signal strength, machine name, etc.  We don't have the file
  1271.     // offsets for these yet, so we must get them for each item.  Because
  1272.     // they're optional, errors are ignored after this point.
  1273.     //
  1274.     if (DEBUG)    cerr << "ReadABI: 5" << endl;
  1275.     // Get Signal Strength Offset
  1276.     if (getIndexEntryLW(fp,indexO,SIGNALtag,1,5,&signalO))
  1277.         {
  1278.         int2 C,A,G,T;
  1279.         int2* base[4];
  1280.         base[0] = &C;
  1281.         base[1] = &A;
  1282.         base[2] = &G;
  1283.         base[3] = &T;
  1284.         if (fseek(fp, signalO, SEEK_SET) >= 0 &&
  1285.             readABIInt2(fp, base[baseIndex((char)(fwo_>>24&255))]) &&
  1286.             readABIInt2(fp, base[baseIndex((char)(fwo_>>16&255))]) &&
  1287.             readABIInt2(fp, base[baseIndex((char)(fwo_>>8&255))]) &&
  1288.             readABIInt2(fp, base[baseIndex((char)(fwo_&255))]))
  1289.             {
  1290.             sprintf(line,"avg_signal_strength = C:%d A:%d G:%d T:%d\n",C,A,G,T);
  1291.             strcat(comments,line);
  1292.             }
  1293.         } // signal strength
  1294.  
  1295.     // the following returns absurd results (ie. on order of 10E10)
  1296.     #if (FALSE)
  1297.     // Get the spacing.. it's a float but don't worry yet
  1298.     int4    spacing;
  1299.     if (getIndexEntryLW(fp,indexO,SPACINGtag,1,5,&spacing))
  1300.         {
  1301.         sprintf(line,"avg_spacing = %lf\n",spacing);
  1302.         strcat(comments,line);
  1303.         }
  1304.     #endif
  1305.  
  1306.     // Get primer position
  1307.     if (getIndexEntryLW(fp,indexO,PPOStag,1,5,&ppos))
  1308.         {
  1309.         // ppos stored in MBShort of pointer
  1310.         primerPosition = ppos>>16;
  1311.  
  1312.         sprintf(line,"primer_position = %d\n",primerPosition);
  1313.         strcat(comments,line);
  1314.         }
  1315.  
  1316.     // machine name
  1317.     if ((getIndexEntryLW(fp,indexO,MCHNtag,1,5,&MCHN_O)) &&
  1318.         (fseek(fp, MCHN_O, SEEK_SET) >= 0))
  1319.         {
  1320.         unsigned char l;
  1321.         char buffer[256];                    // Pascal-type string; len=buffer[0]
  1322.         fread(&l,sizeof(char),1,fp);
  1323.         fread(buffer,l,1,fp);
  1324.         sprintf(line,"machine_name = %.*s\n",l,buffer);
  1325.         strcat(comments,line);
  1326.         }
  1327.  
  1328.     // Dye Primer Offset
  1329.     if ((getIndexEntryLW(fp,indexO,PDMFtag,1,5,&PDMF_O)) &&
  1330.         (fseek(fp, PDMF_O, SEEK_SET) >= 0))
  1331.         {
  1332.         unsigned char l;
  1333.         char buffer[256];                    // Pascal-type string; len=buffer[0]
  1334.         fread(&l,sizeof(char),1,fp);
  1335.         fread(buffer,l,1,fp);
  1336.         sprintf(line,"dye_primer = %.*s\n",l,buffer);
  1337.         strcat(comments,line);
  1338.         }
  1339.  
  1340.     // sample name
  1341.     if ((getIndexEntryLW(fp,indexO,SMPLtag,1,5,&SMPL_O)) &&
  1342.         (fseek(fp, SMPL_O, SEEK_SET) >= 0))
  1343.         {
  1344.         unsigned char l;
  1345.         char buffer[256];                    // Pascal-type string; len=buffer[0]
  1346.         fread(&l,sizeof(char),1,fp);
  1347.         fread(buffer,l,1,fp);
  1348.         sprintf(line,"sample_name = %.*s\n",l,buffer);
  1349.         strcat(comments,line);
  1350.         }
  1351.  
  1352.     if (DEBUG)    cerr << "ReadABI: Exiting normally" << endl;
  1353.     return noError;
  1354.     }
  1355.  
  1356. //    ============================================================================
  1357. //    WriteABI [private]
  1358. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  1359. //    Write a ABI file; called exclusively by Write(FILE*,format_t).
  1360. //    ========================================|===================================
  1361. CTraceFile::error_t
  1362. CTraceFile::WriteABI(FILE* fp)
  1363.     {
  1364.     fp = NULL;
  1365.     return fmtNotSuppError;
  1366.     }
  1367.  
  1368.  
  1369. //    ============================================================================
  1370. //       #    #       #######           ###         # #######
  1371. //      # #   #       #                  #         #  #     #
  1372. //     #   #  #       #                  #        #   #     #
  1373. //    #     # #       #####              #       #    #     #
  1374. //    ####### #       #                  #      #     #     #
  1375. //    #     # #       #                  #     #      #     #
  1376. //    #     # ####### #                 ###   #       #######
  1377. //
  1378. //    ALF format I/O methods
  1379. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  1380. //    
  1381. //    ========================================|===================================
  1382.  
  1383. //    ============================================================================
  1384. //    ReadALF [private]
  1385. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  1386. //    Read a ALF file; called exclusively by Read(FILE*,format_t).
  1387. //    ========================================|===================================
  1388. CTraceFile::error_t
  1389. CTraceFile::ReadALF(FILE* fp)
  1390.     {
  1391.     fp = NULL;
  1392.     return fmtNotSuppError;
  1393.     }
  1394.  
  1395. //    ============================================================================
  1396. //    WriteALF [private]
  1397. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  1398. //    Write a ALF file; called exclusively by Write(FILE*,format_t).
  1399. //    ========================================|===================================
  1400. CTraceFile::error_t
  1401. CTraceFile::WriteALF(FILE* fp)
  1402.     {
  1403.     fp = NULL;
  1404.     return fmtNotSuppError;
  1405.     }
  1406.